clear;
N = 8;
h = 1 / N;
nu = 10;
xh = linspace(0, 1, N + 1);
fh = source(xh) * h^2;
uh = zeros(size(xh));
uh(1) = solution(0);
uh(end) = solution(1);
x2h = linspace(0, 1, N / 2 + 1);
u2h = zeros(size(x2h));
for i = 1:nu
    uh = jacobi(uh, fh);
end
%error = norm(uh - solution(xh))
%plot(xh, uh);
%hold on
%plot(xh, uh - solution(xh), 'r');
rh = fh - Aop(uh);
r2h = fwr(rh)
%plot (xh, rh, 'b*-');
%hold on;
%plot (x2h, r2h, 'ro-');
e2h = zeros(size(r2h));
for i = 1:4000
    e2h = jacobi(e2h, r2h);
end
eh = intp(e2h);
plot (xh, uh, 'b-');
hold on;
uh = uh + eh;
for i = 1:nu
    uh = jacobi(uh, fh);
end
plot (xh, uh, 'ro-');
